{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "</center>\n",
    "Автор материала: программист-исследователь Mail.ru Group, старший преподаватель Факультета Компьютерных Наук ВШЭ Юрий Кашницкий. Сделано на основе курса \"Построение выводов по данным\" специализации \"Машинное обучение и анализ данных\" Яндекса и МФТИ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Тема 9. Анализ временных рядов в Python</center>\n",
    "## <center>Часть 2. Смерти из-за несчастного случая в США</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Известно ежемесячное число смертей в результате случайного случая в США с января 1973 по декабрь 1978, необходимо построить прогноз на следующие 2 года."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = 12, 10\n",
    "from itertools import product\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from scipy import stats\n",
    "\n",
    "\n",
    "def invboxcox(y, lmbda):\n",
    "    if lmbda == 0:\n",
    "        return np.exp(y)\n",
    "    else:\n",
    "        return np.exp(np.log(lmbda * y + 1) / lmbda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deaths = pd.read_csv(\n",
    "    \"../../data/accidental-deaths-in-usa-monthly.csv\",\n",
    "    index_col=[\"Month\"],\n",
    "    parse_dates=[\"Month\"],\n",
    ")\n",
    "deaths.rename(\n",
    "    columns={\"Accidental deaths in USA: monthly, 1973 ? 1978\": \"num_deaths\"},\n",
    "    inplace=True,\n",
    ")\n",
    "deaths[\"num_deaths\"].plot()\n",
    "plt.ylabel(\"Accidental deaths\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проверка стационарности и STL-декомпозиция ряда:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sm.tsa.seasonal_decompose(deaths[\"num_deaths\"]).plot()\n",
    "print(\n",
    "    \"Критерий Дики-Фуллера: p=%f\" % sm.tsa.stattools.adfuller(deaths[\"num_deaths\"])[1]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Стационарность"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Критерий Дики-Фуллера не отвергает гипотезу нестационарности, но небольшой тренд остался. Попробуем сезонное дифференцирование; сделаем на продифференцированном ряде STL-декомпозицию и проверим стационарность:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deaths[\"num_deaths_diff\"] = deaths[\"num_deaths\"] - deaths[\"num_deaths\"].shift(12)\n",
    "sm.tsa.seasonal_decompose(deaths[\"num_deaths_diff\"][12:]).plot()\n",
    "print(\n",
    "    \"Критерий Дики-Фуллера: p=%f\"\n",
    "    % sm.tsa.stattools.adfuller(deaths[\"num_deaths_diff\"][12:])[1]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Критерий Дики-Фуллера отвергает гипотезу нестационарности, но полностью избавиться от тренда не удалось. Попробуем добавить ещё обычное дифференцирование:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deaths[\"num_deaths_diff2\"] = deaths[\"num_deaths_diff\"] - deaths[\n",
    "    \"num_deaths_diff\"\n",
    "].shift(1)\n",
    "sm.tsa.seasonal_decompose(deaths[\"num_deaths_diff2\"][13:]).plot()\n",
    "print(\n",
    "    \"Критерий Дики-Фуллера: p=%f\"\n",
    "    % sm.tsa.stattools.adfuller(deaths[\"num_deaths_diff2\"][13:])[1]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Гипотеза нестационарности уверенно отвергается, и визуально ряд выглядит лучше — тренда больше нет. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Подбор модели"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на ACF и PACF полученного ряда:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.subplot(211)\n",
    "sm.graphics.tsa.plot_acf(\n",
    "    deaths[\"num_deaths_diff2\"][13:].values.squeeze(), lags=58, ax=ax\n",
    ")\n",
    "ax = plt.subplot(212)\n",
    "sm.graphics.tsa.plot_pacf(\n",
    "    deaths[\"num_deaths_diff2\"][13:].values.squeeze(), lags=58, ax=ax\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Начальные приближения: Q=2, q=1, P=2, p=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ps = range(0, 3)\n",
    "d = 1\n",
    "qs = range(0, 1)\n",
    "Ps = range(0, 3)\n",
    "D = 1\n",
    "Qs = range(0, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters = product(ps, qs, Ps, Qs)\n",
    "parameters_list = list(parameters)\n",
    "len(parameters_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "results = []\n",
    "best_aic = float(\"inf\")\n",
    "\n",
    "\n",
    "for param in parameters_list:\n",
    "    # try except нужен, потому что на некоторых наборах параметров модель не обучается\n",
    "    try:\n",
    "        model = sm.tsa.statespace.SARIMAX(\n",
    "            deaths[\"num_deaths\"],\n",
    "            order=(param[0], d, param[1]),\n",
    "            seasonal_order=(param[2], D, param[3], 12),\n",
    "        ).fit(disp=-1)\n",
    "    # выводим параметры, на которых модель не обучается и переходим к следующему набору\n",
    "    except ValueError:\n",
    "        print(\"wrong parameters:\", param)\n",
    "        continue\n",
    "    aic = model.aic\n",
    "    # сохраняем лучшую модель, aic, параметры\n",
    "    if aic < best_aic:\n",
    "        best_model = model\n",
    "        best_aic = aic\n",
    "        best_param = param\n",
    "    results.append([param, model.aic])\n",
    "\n",
    "warnings.filterwarnings(\"default\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Если в предыдущей ячейке возникает ошибка, убедитесь, что обновили statsmodels до версии не меньше 0.8.0rc1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_table = pd.DataFrame(results)\n",
    "result_table.columns = [\"parameters\", \"aic\"]\n",
    "print(result_table.sort_values(by=\"aic\", ascending=True).head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Лучшая модель:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(best_model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Её остатки:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(211)\n",
    "best_model.resid[13:].plot()\n",
    "plt.ylabel(u\"Residuals\")\n",
    "\n",
    "ax = plt.subplot(212)\n",
    "sm.graphics.tsa.plot_acf(best_model.resid[13:].values.squeeze(), lags=48, ax=ax)\n",
    "\n",
    "print(\"Критерий Стьюдента: p=%f\" % stats.ttest_1samp(best_model.resid[13:], 0)[1])\n",
    "print(\n",
    "    \"Критерий Дики-Фуллера: p=%f\" % sm.tsa.stattools.adfuller(best_model.resid[13:])[1]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Остатки несмещены (подтверждается критерием Стьюдента), стационарны (подтверждается критерием Дики-Фуллера и визуально), неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой).\n",
    "Посмотрим, насколько хорошо модель описывает данные:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deaths[\"model\"] = best_model.fittedvalues\n",
    "deaths[\"num_deaths\"].plot()\n",
    "deaths[\"model\"][13:].plot(color=\"r\")\n",
    "plt.ylabel(\"Accidental deaths\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Прогноз"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dateutil.relativedelta import relativedelta\n",
    "\n",
    "deaths2 = deaths[[\"num_deaths\"]]\n",
    "date_list = [\n",
    "    pd.datetime.strptime(\"1979-01-01\", \"%Y-%m-%d\") + relativedelta(months=x)\n",
    "    for x in range(0, 24)\n",
    "]\n",
    "future = pd.DataFrame(index=date_list, columns=deaths2.columns)\n",
    "deaths2 = pd.concat([deaths2, future])\n",
    "deaths2[\"forecast\"] = best_model.predict(start=72, end=100)\n",
    "\n",
    "deaths2[\"num_deaths\"].plot(color=\"b\")\n",
    "deaths2[\"forecast\"].plot(color=\"r\")\n",
    "plt.ylabel(\"Accidental deaths\");"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
